library(dplyr)
library(knitr)
library(maptools)
library(TreeSegmentation)
library(ggplot2)
library(rgl)
library(clue)
library(lidR)
knit_hooks$set(webgl = hook_webgl)

opts_chunk$set(warning=F,message=F)

#set color ramp for treeID
col = pastel.colors(200)

Load in ground-truth

shps<-list.files("/Users/ben/Dropbox/Weecology/ECODSEdataset/Task1/ITC/",pattern=".shp",full.names = T)
itcs<-lapply(shps,readShapePoly)
itcs[[1]]
## class       : SpatialPolygonsDataFrame 
## features    : 9 
## extent      : 402390.1, 402416, 3286283, 3286317  (xmin, xmax, ymin, ymax)
## coord. ref. : NA 
## variables   : 3
## names       : crown_id, confidence,  Plot_ID 
## min values  :       31,       High, OSBS_001 
## max values  :      565,       High, OSBS_001
names(itcs)<-sapply(itcs,function(x){
  id<-unique(x$Plot_ID)
  })
print(names(itcs))
##  [1] "OSBS_001" "OSBS_003" "OSBS_006" "OSBS_007" "OSBS_008" "OSBS_009"
##  [7] "OSBS_010" "OSBS_011" "OSBS_014" "OSBS_015" "OSBS_016" "OSBS_017"
## [13] "OSBS_018" "OSBS_019" "OSBS_025" "OSBS_026" "OSBS_029" "OSBS_030"
## [19] "OSBS_032" "OSBS_033" "OSBS_034" "OSBS_035" "OSBS_036" "OSBS_037"
## [25] "OSBS_038" "OSBS_042" "OSBS_043" "OSBS_044" "OSBS_048" "OSBS_051"
plot(ground_truth<-itcs[[1]])

itcs[[1]]
## class       : SpatialPolygonsDataFrame 
## features    : 9 
## extent      : 402390.1, 402416, 3286283, 3286317  (xmin, xmax, ymin, ymax)
## coord. ref. : NA 
## variables   : 3
## names       : crown_id, confidence,  Plot_ID 
## min values  :       31,       High, OSBS_001 
## max values  :      565,       High, OSBS_001
#proj_itc <-spTransform(itc, CRS("+init=epsg:32617"))

Find corresponding tile

fname<-get_tile_filname(itcs[[1]])
tile<-readLAS(paste("/Users/ben/Dropbox/Weecology/NEON/cropped_",fname,sep=""))
plot(tile)

You must enable Javascript to view this page properly.

Silva 2016 Crown Segmentation Model

silva<-silva2016(tile=tile,extra=T)
## [1] "Computing Ground Model"
## [1] "Computing Canopy Model"
## [1] "Clustering Trees"
##    user  system elapsed 
##   2.432   0.030   2.478 
## [1] "Creating tree polygons"

You must enable Javascript to view this page properly.

Overlay ground truth and predictions

plot(ground_truth,col='red')
plot(silva$silva_convex,add=T)

Classified lidar cloud versus predicted polygons

plot(silva$silva_tile,color="treeID",col=col,size=2)
show2d(face='z+',z=0,{
  plot(silva$silva_convex,col=rgb(255, 0, 0, 30, maxColorValue=255) )
})

Classified lidar cloud versus ground_truth

plot(silva$silva_tile,color="treeID",col=col)
show2d(face='z+',z=0,{
  plot(ground_truth,col=rgb(255, 0, 0, 30, maxColorValue=255) )
})

You must enable Javascript to view this page properly.

Assign Trees

Each tree is assigned based on the maximum overlap. Pairwise membership is done using a Hungarian Algorithm. See clue::solve_LSAP.

assignment<-assign_trees(ground_truth,prediction=silva$silva_convex)

Intersection over union

#loop through assignments and get jaccard statistic for each assignment pair
statdf<-calc_jaccard(assignment=assignment,ground_truth = ground_truth,prediction=silva$silva_convex)
qplot(statdf$IoU)

mean(statdf$IoU)
## [1] 8.476192
median(statdf$IoU)
## [1] 0.1739069